!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Semidiscretized ODE for the compressible Navier-Stokes DG system
!!
!! \n
!!
!! This module represents the semidiscretized problem as an ODE, with
!! the layout required by the time integrators of \c
!! mod_time_integrators.
!<
module mod_dgcomp_ode

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_time_integrators, only: &
   mod_time_integrators_initialized, &
   c_ode

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_bcs, only: &
   mod_bcs_initialized, &
   t_bcs

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_dgcomp_rhs, only: &
   mod_dgcomp_rhs_initialized, &
   t_bcs_error, &
   dgcomp_tens

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_dgcomp_ode_constructor, &
   mod_dgcomp_ode_destructor,  &
   mod_dgcomp_ode_initialized, &
   t_dgcomp_ode,          &
   new_dgcomp_ode, clear

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 
 !> ODE type
 !!
 !! The state fields are allocatable components, to allow different
 !! variables to represent diffent states. The fields that represent
 !! the underlying discretization, on the other hand, are pointers,
 !! since all the variables have to point to the same configuration.
 type, extends(c_ode) :: t_dgcomp_ode
  type(t_grid), private, pointer :: grid => null()
  type(t_base), private, pointer :: base => null()
  type(t_bcs ), private, pointer :: bcs  => null()
  logical, private :: viscous_flow
  real(wp), allocatable :: u(:,:,:)
  type(t_bcs_error) :: bcerr
 contains
  procedure, pass(tnd) :: rhs => dgcomp_ode_rhs
  procedure, pass(z)   :: add => dgcomp_ode_add
  procedure, pass(y)   :: stimes => dgcomp_ode_stimes
  ! no need to overload solve since we use explicit solvers
 end type t_dgcomp_ode

 ! private members

 interface clear
   module procedure clear_dgcomp_ode
 end interface

! Module variables

 ! public members
 logical, protected ::               &
   mod_dgcomp_ode_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_dgcomp_ode'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_dgcomp_ode_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
(mod_time_integrators_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
         (mod_bcs_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) .or. &
  (mod_dgcomp_rhs_initialized.eqv..false.) ) then 
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_dgcomp_ode_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_dgcomp_ode_initialized = .true.
 end subroutine mod_dgcomp_ode_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_dgcomp_ode_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_dgcomp_ode_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_dgcomp_ode_initialized = .false.
 end subroutine mod_dgcomp_ode_destructor

!-----------------------------------------------------------------------
 
 subroutine new_dgcomp_ode(obj,grid,base,bcs,viscous_flow,nvar)
  type(t_grid), intent(in), target :: grid
  type(t_base), intent(in), target :: base
  type(t_bcs ), intent(in), target :: bcs
  logical, intent(in) :: viscous_flow
  integer, intent(in) :: nvar
  type(t_dgcomp_ode), intent(out) :: obj

   obj%grid => grid
   obj%base => base
   obj%bcs  => bcs
   obj%viscous_flow = viscous_flow
   allocate( obj%u(nvar,base%pk,grid%ne) )
   ! leaving the default initialization of bcerr

 end subroutine new_dgcomp_ode
 
!-----------------------------------------------------------------------
 
 pure subroutine clear_dgcomp_ode(obj)
  type(t_dgcomp_ode), intent(out) :: obj

   nullify( obj%grid )
   nullify( obj%base )
   nullify( obj%bcs  )
   ! allocatable components are implicitly deallocated
 
 end subroutine clear_dgcomp_ode
 
!-----------------------------------------------------------------------

 subroutine dgcomp_ode_rhs(tnd,t,uuu)
  real(wp),            intent(in)    :: t   !< time level
  class(c_ode),        intent(in)    :: uuu !< present state
  class(t_dgcomp_ode), intent(inout) :: tnd !< tendency

   select type(uuu); type is(t_dgcomp_ode)

   call dgcomp_tens( tnd%u , uuu%grid , uuu%base , uuu%bcs , uuu%u , &
                     tnd%bcerr , uuu%viscous_flow )

   end select
 end subroutine dgcomp_ode_rhs

!-----------------------------------------------------------------------

 pure subroutine dgcomp_ode_add(z,x,y)
  class(c_ode),        intent(in)    :: x, y
  class(t_dgcomp_ode), intent(inout) :: z

   select type(x); type is(t_dgcomp_ode)
     select type(y); type is(t_dgcomp_ode)

   z%u = x%u + y%u

   ! propagate bcs errors
   if(x%bcerr%lerr) then
     z%bcerr%lerr    = .true.
     z%bcerr%message = x%bcerr%message
   elseif(y%bcerr%lerr) then
     z%bcerr%lerr    = .true.
     z%bcerr%message = y%bcerr%message
   else
     z%bcerr%lerr    = .false.
   endif

     end select
   end select
 end subroutine dgcomp_ode_add

!-----------------------------------------------------------------------

 pure subroutine dgcomp_ode_stimes(y,r,x)
  real(wp),            intent(in)    :: r
  class(c_ode),        intent(in)    :: x
  class(t_dgcomp_ode), intent(inout) :: y

   select type(x); type is(t_dgcomp_ode)

   y%u = r*x%u

   y%bcerr = x%bcerr

   end select
 end subroutine dgcomp_ode_stimes

!-----------------------------------------------------------------------

end module mod_dgcomp_ode

